Using force covariance to derive effective stochastic interactions in dissipative particle dynamics 
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There exist methods for determining effective conservative interactions in coarse grained particle based meso- 
scopic simulations. The resulting models can be used to capture thermal equilibrium behavior, but in the model 
system we study do not correctly represent transport properties. In this article we suggest the use of force covari- 
ance to determine the full functional form of dissipative and stochastic interactions. We show that a combination 
of the radial distribution function and a force covariance function can be used to determine all interactions in 
dissipative particle dynamics. Furthermore we use the method to test if the effective interactions in dissipa- 
tive particle dynamics (DPD) can be adjusted to produce a force covariance consistent with a projection of a 
microscopic Lennard- Jones simulation. The results indicate that the DPD ansatz may not be consistent with 
the underlying microscopic dynamics. We discuss how this result relates to theoretical studies reported in the 
literature. 
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I. INTRODUCTION 

It is the long-standing aim of molecular simulations to elu- 
cidate mechanisms that cannot be directly observed in ex- 
periments, or understood in terms of more abstract models. 
Though extremely successful in many areas, when applied 
to mesoscopic systems, such as membranes or complex flu- 
ids, one often finds that the relevant time and length scales 
on which important mechanisms take place are beyond the 
reach of direct, detailed, simulations. As a consequence, sev- 
eral coarse-grained methods have been developed to allow for 
a larger time span to be simulated. Lattice gases (jl corre- 
sponds to dividing the system into an array of sub-systems, 
each a thermodynamic system on its own with a local tem- 
perature, pressure, particle density and velocity distribution. 
Other coarse-graining procedures have explicit particles with 
pair-wise interactions; well-known examples are united atoms 
yfl, smoothed particle hydrodynamics (SPH) y[| and dissipa- 
tive particle dynamics (DPD) 0]. There also exists hybrid 
methods suitable for example in molecular simulations where 
atomistic resolution is needed only in spatially localized do- 
mains, e.g. JHH]. 

Atomic force fields for molecular dynamics (MD), derived 
from potentials defined empirically or theoretically (e.g. from 
quantum-mechanical models), are relatively mature. In con- 
trast, it is much less clear how to choose the effective force 
fields for coarse-grained models, partly because the connec- 
tion between the degrees of freedom in the coarse-grained dy- 
namics and the underlying molecular dynamics differ from 
one coarse-graining procedure to the next. Frequently, sim- 
ple heuristic forces are used; partly because of computational 
ease but also because the detailed forces may not be known 
0]. The magnitude of the forces are then chosen to match 
macroscopic observables of the system, such as the compress- 
ibility. Forrest and Suter |8(| (see also references therein) cal- 
culate an effective force on the particles, interacting according 
to the Lennard- Jones potential, which correspond to the aver- 
age effect of the true forces during a time interval. Because the 
time-averaged force is effectively an average over rapid fluc- 
tuations of close particles, the effective potentials are much 



softer than the Lennard-Jones potential at small distances. An 
alternative approach, which we will use later in this paper, 
is to use the fact that for particle systems with central forces 
there is a one-to-one relation between the radial distribution 
of particles at thermal equilibrium and the pair-wise potential 
& 

When coarse-graining a molecular system, the effective in- 
teractions in the resulting system can either be determinis- 
tic (smooth, due to averaging of fast degrees of freedom), in 
which case just matching the equilibrium properties of the sys- 
tem gives the correct dynamics, or the fast degrees of freedom 
act as a noisy driving force. Which of these two scenarios 
that best describe the system at hand depends on the exchange 
of energy between the simulated system and the degrees lost 
in the coarse-graining procedure. If a substantial amount of 
energy is exchanged then the motion of the coarse-grained 
particles will not be smooth or deterministic. In this case, 
it is not sufficient to capture the conservative forces but we 
must also introduce dissipative and stochastic forces. These 
forces are included in the SPH and DPD models. The mod- 
els are quite similar, and in this article we choose to focus on 
DPD. As mentioned in the previous paragraph the smooth, or 
conservative, part of the interaction can be determined from 
the radial distribution function in in thermal equilibrium @]. 
In other words, the radial distribution at equilibrium is not 
affected by the noise term in the interaction (provided an ap- 
propriate deterministic dissipative force is added, ensuring the 
equilibrium temperature). Clearly it follows that equilibrium 
properties that are determined by the radial distribution are 
also not affected. Examples are compressibility and other ob- 
servables defining the equation of state. Other characteristics 
such as the diffusion coefficient and the viscosity do depend 
on the stochastic interactions. To capture these properties it is 
central to choose the stochastic interaction, including its radial 
dependence, as correctly as possible. 

To accurately describe the dissipative and stochastic part of 
the dynamics we must introduce a new observable, comple- 
menting the radial distribution function. A candidate could be 
the autocorrelation of the velocity, which is directly related to 
the stochastic driving on a single particle as well as the dif- 
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fusion coefficient through Green-Kubo relations 
a particle-based mesoscopic model such as DPD, the stochas- 
tic interaction is represented as a force between the particles. 
Clearly this interaction must have a radial dependence (typi- 
cally with short range). A useful observable for estimating the 
stochastic interaction must therefore be able to resolve the ra- 
dial dependence. In this paper we suggest the use of the force 
covariance function as a candidate for such an observable, and 
use DPD as a test case for the method. 



II. DISSIPATIVE PARTICLE DYNAMICS 

DPD was introduced in 1992 by Hoogerbrugge and Koel- 
man |0] as a simulation technique for hydrodynamic phe- 
nomena. The method has received much theoretical atten- 
tion iH [3 [M [3 [Hi that provides support for this kind 
of model and has been established as a standard method for 
mesoscopic simulation. Amongst other things DPD has been 
used to study complex fluids Q, spontaneous self-assembly 
of amphophilic molecules into bilayered membranes fl7il . 
vesicles fl8l fl9tl . and hydrodynamics l20ll . In its standard 
form, DPD is a particle model with pairwise interactions, 
quite similar to molecular dynamics, but with a dissipative and 
stochastic contribution to the interactions between the parti- 
cles. 

In its simplest form, the equations of motion for a DPD 
model, with mesoscopic particles positioned at r;, with ve- 
locities Vj and momenta can be written as a system of 
Langevin equations 

= Vi, 

»>< E[ F £ ' v " ' p «]. CD 

where F^ , F° and are the conservative, dissipative and 
stochastic forces between particles i and j. Both the conser- 
vative and non-conservative interactions in DPD are modeled 
by central forces obeying Newton's third law, ensuring that 
(angular) momentum is conserved @J. The dissipative and 
stochastic forces are 

F°. = -uj d ( rij ) e 4J • ( Vi - Vj ) eij , (2) 

F y' = wS ( r «) Cij e 'ji (3) 

where is the distance between particles i and j, and e,j 
is the unit vector pointing from j to i. The scalar functions 
up (r^) and ui s (rij) describe friction and noise, respectively. 
Qij is interpreted as a symmetric Gaussian white noise term 
with mean zero and covariance 

(dj{t)Q'j'(t')) = (*«'%' + Sa>S JV )6{t - 0, (4) 

where <5y and S(t) are the Kronecker and Dirac delta func- 
tions. Assuming that the equilibrium distribution of a DPD 
system is given by the canonical ensemble, the fluctuation- 
dissipation theorem leads to a relation between the dissipative 
and stochastic parts 11211 : 

u D (r) = (2k B T)- 1 [co s (r)} 2 . (5) 



For simplicity, we drop the superscript and write ui(r) = 
u) s (r). Eqs. (I])-© together establishes the general form of 
the DPD dynamics. Both the conservative force F^ , or equiv- 
alently the corresponding scalar potential, and the scalar func- 
tion uj(r) depend on the particular system of interest and need 
to be determined to obtain the correct DPD model. In prac- 
tice, this is the difficult part of DPD, and also the rationale 
behind the heuristic approach in deciding the interactions. As 
an example, the common practice for fluid like systems is to 
apply linear functions with a cutoff radius r c ; 

Vfj = (1 - nj /r c ) a,, Xij ey = / ' (rtj ) e,, , (6) 

u(rij) = (1 - rij/r c )axij, (7) 

where By is the strength of the conservative force between 
particles i and j, a is the amplitude of the noise, and Xij is 
one for < r c and zero elsewhere. 

Despite its popularity and theoretical support, it is un- 
clear how DPD should be interpreted as a coarse-grained 
model [21]. One point of view, and the one we will elab- 
orate on in this article, is to consider DPD as a systematic 
coarse-graining of an underlying atomistic system. If the DPD 
method could be shown to have a firm microscopic founda- 
tion, that would greatly impact our ability to couple DPD to 
actual physical systems. Several authors, e.g. j22il23ll24ll25ll . 
have established bottom-up connections between the micro- 
and mesoscale and obtained mesoscopic dynamics resembling 
DPD. The resulting methods differ from DPD as they incor- 
porate the geometry of the system in the equations, implying 
forces that are not central or pairwise, while DPD is a model 
with only pairwise interactions. It is clear that the validity 
of DPD as a coarse-grained model, or how well DPD repre- 
sents an underlying microscopic system, has not been fully 
resolved. 

To obtain a well-defined bottom-up scheme, the dynamics 
of the coarse-grained DPD particles must be defined through 
a projection of the microscopic trajectories. The problem is to 
find a closed representation of the system at the coarse-grained 
level, i.e. to determine all interactions in the DPD model. In 
this article we investigate a method of estimating the DPD in- 
teractions using measurements on the coarse-grained level of 
a simulation. By applying the method to a typically assumed 
projection of a microscopic system, we clarify some important 
aspects of DPD as a systematically coarse-grained model. 

The DPD technique has its theoretical foundations in Mori- 
Zwanzig theory on projection operators rf26l I27I I2H, l29ll . In 
short, the theory states that given a microscale dynamics, 
a lower dimensional representation can be formally attained 
through a projection of the phase space, where fast degrees 
of freedom are treated as Markovian (white) noise 12611. This 
framework can be applied to molecular dynamics l30l l3lll . 
Naturally, how faithfully the coarse-grained model will repre- 
sent the underlying dynamics depends on the choice of pro- 
jection. The DPD method assumes a projection resulting in a 
mesoscopic model characterized as a particle based Langevin 
dynamics with pairwise and negated central forces. The inter- 
nal degrees of freedom in the mesoscopic particles give rise to 
dissipation and noise, which is captured by non-conservative 
pairwise interactions. As a consequence, sufficiently close to 
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equilibrium, one obtains the classical result of the asymptotic 
t~ d / 2 decay of the velocity auto-correlation (d is the dimen- 
sionality of the system) 0211 . In addition, the interactions give 
rise to hydrodynamic modes in the fluid lf32T[33ll . which lead 
to the Navier-Stokes equations on the macroscopic level 



III. ESTIMATING THE EFFECTIVE FORCES 

Given the DPD ansatz for the effective equations of motion, 
the question is: how does one find the conservative and dissi- 
pative forces F c (r) and u>(r)l In this section we present the 
theoretical motivations for our method, and apply it to DPD 
simulations to test the accuracy of the method on a case where 
we know the ansatz to be true. In section [IV] we apply the 
method to a coarse-graining of a system of particles interact- 
ing via the Lennard- Jones potential, in order to see how the 
method fares on a classical molecular dynamics system. 



A. The conservative force term 

The original motivation §j\\ for a repulsive conservative 
force was a measurement of the effective potential for the 
interaction between particles in a Lennard- Jones fluid @|. 
More generally applicable methods for estimating the conser- 
vative interactions are based on the radial distribution function 
(RDF) in thermal equilibrium il IMS [Mill III . In these 
reports, the estimate of the conservative force is obtained us- 
ing a result by Henderson U, stating that the difference be- 
tween two pairwise potentials that give rise to the same radial 
distribution function must be a constant gauge shift, and hence 
of no physical significance. The importance of this theorem 
lies in the one-to-one correspondence between potential and 
radial distribution function. 

The conservative interactions are determined by the RDF 
alone, which in turn is determined by the thermal equilibrium 
of the system. As long as the fluctuation-dissipation theorem 
holds, the thermal equilibrium is independent of the specific 
form of the dissipative and random interactions fl2ll . and it 
follows that we can estimate the conservative forces from a 
given RDF independently of the stochastic forces. Here we 
use the inverse Monte Carlo method of Lyubartsev and Laak- 
sonen l4lll . which starts from a Boltzmann ansatz of the po- 
tential and then, through iteration, finds a potential giving rise 
to the desired RDF. In what follows we briefly describe the 
method, following fiHl . 

The connection between the RDF and the potential can be 
found from the Hamiltonian of the system. Consider a system 
of particles with pairwise interactions. It can be discretized as 



radial distribution function, g(r), by the relation 



(8) 



which corresponds to using a stepwise constant potential, <& Q . 
S a denotes the number of particle pairs separated by a dis- 
tance in the range r a to r a+ i, where tq = 0, ?'i = dr and 
r a = a ■ dr. The average of S a is directly connected to the 



(Sa) 



v a 



N(N -l)/2 L 3 



(9) 



where N is the number of particles, L 3 the volume of the sim- 
ulation box and V a the volume of the spherical shell between 
radii r a and r a+ \. Using a Monte Carlo (MC) approach, the 
system may be simulated with a start potential ■ It is com- 
mon practice to choose this to be the potential of mean force, 



$L 0) =-k B T\ng{r a ). 



(10) 



The correspondance between potential and RDF is an equilib- 
rium result and hence only valid for fixed temperatures and 
densities. These quantities must therefore be the same in the 
MC simulation as they were in the original simulation from 
which the RDF was obtained. 

Simulating with the trial potential, , produces an (S'i ' 1 ) 
which may differ from the correct value, S**. The difference, 
A(5" Q ) (0) = (S Q ) (0) - S*, is used to find a new trial potential 
by solving for A$ in the linear equation system 



(11) 



with 



given by |41J 



d(S a ) (S a S 7 ) - (5„)<S 7 ) 



k B T 



(12) 



The next guess for potential is then $ < - 1 ' = $ < -°- ) — A$. This 
procedure is repeated until $ has converged to a potential that 
reproduces the original RDF. 

The potential of mean force is usually a good first approxi- 
mation to the final potential, and convergence to a unique po- 
tential normally takes less than ten updates in the Monte Carlo 
simulations. This is especially true for the soft coarse-grained 
potentials we get from considering effective interactions be- 
tween clusters of particles. For instances that nevertheless 
require special care, problems with convergence for the po- 
tential over successive MC simulations can generally be over- 
come by moving only a fraction in the direction specified by 
Eq. O- 



B. The dissipative force term 

Assuming that the DPD ansatz is valid, the functional form 
of the dissipative term (and through the fluctuation dissipation 
theorem, Eq. (O, the stochastic term) can be isolated through 
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a Kramer-Moyal expansion l42tl of Eq. (Q]): 

-{Sp i -Sp j )/St = {u 2 {r ij ))~ ]T (Fg-F%) 6t 

(uj 2 (rik)uJ 2 (rji) (e lk -v ik ) (e^-v^) (e^-e^)) 



E 



(2fc B T)2 



(^ 2 (r^)(e, fc -v a )(e. tfc -F^)) ^ 



2 ( " 2(rj7)(e ^; v t )(ejrF ' fc)) «+q(^), (B) 



2k B T 



where Spi(t) = Pi(t + St) — Pi(t), All averages in Eq. (lE3b 
are conditioned on that the distance r,j between particles i and 
j equals r. This equation provides a relationship between the 
functional form of the stochastic and dissipative interactions, 
to(r), and the force covariance, Kp, defined as: 



k f = -(Spi-6pj)/5t. 



(14) 



In the DPD simulations we can take St small enough that 
only the leading term of Eq. ( fT3l is significant. If the mi- 
croscopic dynamics is deterministic, like in most molecular 
dynamics, there generally exists a time-scale below which the 
forces are smooth functions of time (this is the time-scale on 
which the molecular dynamics can be integrated). On this 
time-scale, the projected dynamics is also smooth (if the pro- 
jection is smooth) but not autonomous. It follows that for 
small St the force covariance Kp is proportional to St, cor- 
responding to ui(r) = in Eq. (fl3t . 

In fluids, the magnitude of ftp will typically increase with 
increasing St, because on relatively short time scales the par- 
ticles of the fluids oscillate in a cage formed by their clos- 
est neighbours (the Franck-Rabinowitch effect). We consider 
these rapid fluctuations to correspond to fast degrees of free- 
dom in the system. To get an idea of the time-scales in- 
volved, consider water particles in a fluid at normal pressure 
and room temperature. The particles' distance to their closest 
neighbours oscillate around the first peak in the radial distri- 
bution function, at approximately 0.28 nm. The half-width of 
the peak, approximately 0.05 nm, gives an indication of how 
far the molecule travels before experiencing strong repulsive 
forces from other particles. We estimate the typical velocity 
as the root mean square (RMS) velocity 



wrms 



(15) 



At room temperature (25 °C), the RMS velocity is approxi- 
mately 640 m/s. One may argue that the orientation of the 
particle velocities are essentially random, so that the RMS 
difference in velocity is i>rms v2> an d that they collide at half 
the half-width. The time to travel this distance at the typical 
velocity is then approximately 0.03 ps, and we take this as a 
rough approximation to the timescale at which the fast dynam- 
ics occur. It is only at timescales significantly larger than this 
timescale that we can expect to approximate the fast degrees 
of freedom with a spatially structured but Markovian noise as 



in the DPD ansatz. On this time scale, the fluid approaches 
a local thermal equilibrium on the length scale of the coarse- 
grained particles, determined by the local concentration, local 
average velocity and kinetic energy l32tl . 

As a consequence, it is generally not possible to take the 
limit of St — > in the numerically estimated ftp to find w(r). 
Rather, we will assume that there exists a time interval where 
the fast degrees of freedom can be approximated by noise, 
and where ftp is an approximately linear function of time. 
Given two values of St in this interval, St\ and St?,, we can use 
Richardson extrapolation to eliminate the St term in Eq. ( fT3] l 
to obtain an 0(St 2 ) estimate for uj 2 (r): 

2/ \i _ Six (Spi-Spj)^ St 2 (Sp l -Sp j )\ 5tl 

w {rij)Ut °~sf 2 st 2 -st, — sf, st 2 -s tl ■ (16) 

An alternative approach is to do a linear fit with respect to St 
in this region, for each value of r, and from the best fit take 
the intersection with the line St = 0. 



C. Recreating the effective interactions of DPD simulations 

At this point we have established the principles behind 
our method. An important consistency check is to apply the 
method to standard DPD simulations, where the dynamics is 
truly Langevinian. This was done by performing DPD simula- 
tions with different functional forms of both (r) and w(r). 
Using standard DPD units, the simulation region was a peri- 
odic cubic box with side length L = 8.7359, with 3 particles 
per volume unit, giving a total of 2000 particles. From the 
simulations, the RDF and Kp were calculated for 100 r-values 
in the range to 1.75, after which the RDF had converged 
to 1.0. The time-step size used in the simulations was small 
(St = 10~ 3 ) compared to a normal DPD simulation. The rea- 
son for this was to approach the limit of small St, so that the 
terms proportional to 5t could be ignored in Eq. ( fT3l ; uj 2 (r) 
is then given simply by «p. 

In all cases the method accurately recreated the DPD in- 
teractions used in the simulations. Fig. Q] shows the results 
of recreating ui 2 {r) for three different functional forms. The 
conservative potential was also varied (see figure caption for 
details), and plots of recreated potentials from these simula- 
tions are shown in figure [2] 

An example of the situation where we cannot measure Kp 
in the limit <5t — > is shown in Fig. [3] Here two measure- 
ments of Kp from a DPD simulation using time steps of differ- 
ent size (Sti = 0.025 and St 2 = 0.05) deviates clearly from 
w 2 (r). The resulting estimate of ui 2 (r), obtained by Richard- 
son extrapolation of Kp measurements, falls close on the orig- 
inal curve. It should be noted that Kp measurements are ob- 
tained from the same simulation (with time-step 5t = 0.005), 
as an increase in the DPD time step would alter the dynamics 
of the system. For a projected dynamics this is not a problem, 
as it evolves on the microscopic time scale. In Fig. [4] mea- 
surements of ftp from the same simulation is plotted against 
the size of the time difference between measurements, St. As 
predicted by Eq. ( fT~3b the system exhibits a linear behaviour 
for small values of St (in this case St < 0.05). Note, however, 
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FIG. 1: The plots show three different functional forms of L0 2 (r). 
The symbols (o, o and □) show the values found by our method: mea- 
suring ftp of a DPD simulation. The exact functional forms used in 
the simulations are plotted as solid lines. All units are standard DPD 
units. The conservative and dissipative forces are (for r £ [0, 1]): 
(A) F c (r) = 10(1 - r), cu(r) = 5r(l - r) (B) F c (r) = 10(1 - r), 
w(r) = 2(1 - r) (C) F c (r) = 10r(l - r), u(r) = 3(1 - r). For 
r > 1, both functional forms are zero. 
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FIG. 2: The conservative potentials $(r) = J°°dr'F (r') from 
three DPD simulations have been recreated from RDF-data. Cases 
(A) and (B) correspond to a linear DPD-force (i.e. quadratic po- 
tential) with different random force parts. For case (C), a quadratic 
conservative force was used. For details, see caption of figure [T] In 
all three cases, the potential was exactly recreated up to statistical 
accuracy. 
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FIG. 3: The plot shows two measurements of kf obtained from a 
DPD simulation using time steps of different size (o: Sti — 0.025 
and o: Sti = 0.05 ). The resulting estimate of Lo 2 (r), obtained by 
Richardson extrapolation of kf measurements, is shown as □. The 
solid line shows the exact form of co 2 (r) = (3(1 — r)) 2 used in the 
simulation. The conservative force used was F c (r) = 10r(l — r). 
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FIG. 4: The force covariance kf measured in DPD simulation plot- 
ted as a function of St, for different values of r (solid lines). The 
simulation setup is the same as in Fig. [3] It is clearly visible that up 
has a linear region for small St (marked by dotted vertical line), as 
expected from Eq. l |13l >. The dashed lines are linear functions with 
slopes given by the derivatives of kf close to St — 0. 



that how far the linear region extends varies significantly with 
the value of r. 



IV. COARSE-GRAINING OF A LENNARD- JONES FLUID 

We now apply the method to a coarse-grained molecular 
system. Because of its simple form and because it is so well 
understood, we examine the case of a coarse-grained two- 
dimensional Lennard-Jones (LJ) fluid. This is a single-species 
fluid, where the particles interact according to the standard 
pairwise potential 

V(r) = 4e [(r/a L3 )- 12 - (r/a u )- 6 } . (17) 



The parameters are chosen to correspond to bulk water at 
room pressure and temperature: the energy e = 6.739 meV, 
the interaction length olj = 0.31655 nm, and the mass 
m LJ = 2.99 x 10~ 26 kg. 

In DPD, the particles are often understood to be a collec- 
tion of underlying particles, with properties such as mass and 
momentum defined from these. According to this view, we 
follow [22], where the coarse-grained dynamics is expressed 
in terms of a set of N mesoscopic particles. Each particle has 
a position Rfc, a velocity Ufe, and a mass Mf.. The instanta- 
neous momentum of mesoscopic particle k is defined as the 
sum of the momenta of the microscopic particles for which k 
is the nearest mesoscopic particle, and the mass of the particle 
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is defined as the total mass of these underlying particles: 

n 

?:=i 

n 

P k = ^^(rOmiVi, (18) 

i=l 

U fc = Rfc = P fe /Af fe . 

Here n is the total number of microscopic particles, and m,-, 
and Vj represent masses, positions and velocities, respectively, 
of the microscopic particles. £fc(i\) is 1 if mesoscopic particle 
k is closer to microscopic particle i than any other mesoscopic 
particle is, and otherwise. 

Though we use Eqs. ( TT~8b to find the motion of the meso- 
scopic particles, it is nevertheless illuminating to see how the 
effective forces acting on the mesoscopic particles are related 
to the forces acting on the microscopic particles. Formally, we 
calculate the time-derivative of the momentum of the meso- 
scopic particle k in Eq. ( fT8l . Between each passage of a mi- 
croscopic particle from one mesoscopic particle to the next, 
£fc(i"i) is constant (either zero or one). During these time in- 
tervals, the effective force acting on the mesoscopic particle 
is the sum of the forces acting on the microscopic particles 
closest to k: 

n 

M k lJ k = J2z k (r i )f i) (19) 

Suppose microscopic particle i leaves mesoscopic particle k. 
When this happens, the mesoscopic particle experiences an 
impulse 

I = AM fc U fc = -nnYi, (20) 

so that the velocity of mesoscopic particle k changes instanta- 
neously from XJk to 

14 = — - (M fc U fc - miVi) . (21) 

The receiving mesoscopic particle is subject to the opposite 
impulse, —I (formally it is possible to express the force in 
terms of Dirac's delta function). 

Finally, a word of caution: It might seem natural to use 
Eq. (T% alone to define the motion of the mesoscopic par- 
ticles; however, in this dynamics the total momentum in 
the mesoscopic system changes when a microscopic particle 
moves from one mesoscopic particle to another. 

A. Estimated forces 

The conservative interaction was determined from the RDF 
of the mesoscopic particles by the inverse MC method dis- 
cussed earlier. The RDF was measured in LJ simulations with 
1600 particles in a simulation box with side length 12.48 nm, 
temperature 333 K and periodic boundaries. In the coarse- 
grained description, 160 particles were used, resulting in an 
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FIG. 5: A: Lennard- Jones potential used to simulate the microscopic 
particles. B: The effective potential for the coarse-grained system, 
obtained using the inverse MC method. C: The standard quadratic 
DPD potential, scaled to the same magnitude as the estimated poten- 
tial. The main characteristics of DPD, soft-core repulsion and finite 
support, are confirmed by the retrieved potential. 



average of 10 microscopic particles per mesoscopic particle. 
Fig. |5] shows the potential compared with both the LJ poten- 
tial ( fTTI i and the standard DPD potential ©. The retrieved po- 
tential confirms the main characteristics of the standard DPD 
potential, i.e. soft-core repulsion and finite support. 

In Fig. [6] we show «p for the Lennard-Jones system as a 
function of St for different values of r (the inset shows Kp 
as a function of r for different values of St). For small St, 
kf is increasing up to a maximum. We find a range starting 
at St = 0.1 ps where Kp is approximately linear. In prin- 
ciple Kp is also linear for very small values of St, but since 
we know that the fluctuations we want to approximate with 
Markovian noise occurs on time scales < 0.03 ps (c.f sec- 
tion HUB), we reject this region. Hence, the linear region in 
the figure should match the linear region of Eq. ( fT3b . An ex- 
trapolation using data from the linear region gives the shape 
of oj 2 (r) shown in Fig. [7] For comparison we have also in- 
cluded the standard shape of w 2 (r), which can be derived 
from Eq. (0, with the magnitude scaled to fit the obtained 
u> 2 (r) (solid line). We note that the dissipative force derived 
from the mesoscopic particle motion is significantly broader 
than the standard shape. The conservative force is increasing 
only gradually as a pair of mesoscopic particles come within 
the interaction distance (approximately 1.5 nm), the dissipa- 
tive force grows much more rapidly. 



B. Consistency check 

To test if the projected LJ-system can be represented by 
pairwise Langevinian dynamics, we perform a DPD simula- 
tion using the estimated functional forms of the conservative 
and dissipative forces, as shown in Fig. [5] and Fig. [7] Setting 
up the DPD simulation so as to correspond to the projected 
LJ-dynamics, we obtained measurements of Kp for varying St. 
In Fig.[8]the results (symbols) for a selection of r-values are 
plotted together with the corresponding measurements from 



7 




0.1 0.15 

St [ps] 
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FIG. 6: The force covariance kf as a function of St (symbols), for 
three values of r (shown by each curve), kf is approximately linear 
for large enough St (to the right of the dashed line). For each value of 
r, the extrapolation of this region (indicated by solid lines) to St — 
determines the values of Lo 2 (r), c.f. Eq. d 1 3 fc - Inset: kf as a function 
of r, for St = 0.1 ps (top), 0.15 ps, 0.2 ps and 0.25 ps (bottom). 
It is clear from this figure that the terms proportional to St are not 
negligible in Eq. J 1 3 b . 




FIG. 7: The plotted circles show the estimate of uj 2 (r), calculated 
from the force covariance kf using Richardson extrapolation, c.f. 
Fig. [6] The solid line is cj 2 (r), using the DPD form in Eq. l(7}, 
with the same cut off distance as for the conservative potential in 
Fig-El The estimated stochastic interaction differs significantly from 
the function commonly used in DPD studies. 



the LJ-projection (solid lines), and the Richardson extrapola- 
tion of these (dashed lines). The linear region of kf for the 
DPD dynamics lies between St = and St m 0.05 ps. As 
pointed out in section HIlB I this is in the region of the fast dy- 
namics for the underlying system. Clearly, the linear regions 
for k f in the DPD system and in the projected Lennard- Jones 
system do not coincide and therefore we cannot confirm that 
the projection can be formulated in terms of DPD. 

As is seen in section IIIIBI the method we have devel- 
oped works for any system that obeys the DPD ansatz. More 
strictly, it works for any system that evolves on a timescale 
where all terms of order 0(St 2 ) can be neglected in Eq. H3i . 
That allows uj(r) to be estimated either directly from k f (if 
also first order terms of St are negligible), or through Richard- 
son extrapolation of Kp for different values of St. If we only 




FIG. 8: The plot shows measurements of the force covariance kf 
plotted as a function of St for both the LJ-projection (solid lines) and 
a DPD simulation (symbols) using the estimated functional forms. 
Also plotted is the Richardson extrapolation of the LJ-projection 
(dashed lines). The curves are plotted for three values of r (same 
as in Fig. |6). The linear region of kf for the DPD dynamics lies 
between St = and St — 0.05 ps, while the linear region of the 
projected dynamics is in the range St = 0.1 to St = 0.25 ps. kf 
for the two systems do not have a coinciding linear region (for each 
value of r). 



have data available for the system on a longer timescale, there 
may not be a region where Kp is approximately linear. In this 
case it is not possible to use the linearization procedure to ex- 
tract the dissipative force. 

An important remark is that, in general, it cannot be con- 
cluded from measurements of kf alone if the dynamics is 
Markovian or follows the DPD ansatz; a cross check with a 
DPD simulation is necessary. In the case of the projected LJ- 
dynamics, it was reasonable to assume that the linear region 
(see Fig. [6} could be interpreted as the right timescale to con- 
sider for extracting the functional forms. However, kf from 
the DPD simulation turned out to have its linear region on 
a much shorter timescale than assumed for the projected dy- 
namics. 

In the light of this result, there are two explanations for the 
observed deviations from the DPD ansatz, which differ with 
respect to whether the DPD ansatz is correct or not. If we 
first assume that the projected system follows the DPD ansatz, 
our guess of a linear region is not correct, and it follows that 
higher order terms of St will affect the value of «f, rendering 
our method unapplicable for this case. The solution to this 
problem calls for more sophisticated methods to estimate oj(r) 
from Kp- The second possibility is that the projection does not 
produce a dynamics that follows the DPD ansatz. This could 
either simply be a result of our choice of projection, or it could 
point to deeper problems with constructing a coarse-graining 
scheme that leads to the DPD model. 

Flekk0y et al. I22I l23l f24ll have used the same type of pro- 
jection as presented in this article, c.f. Eq. ( fl~8T >. but rather 
than considering the coarse grained entities as spherical parti- 
cles, they consider them as cells on a Voronoi lattice. Within 
each cell, the fluid is assumed to correspond to an ideal fluid 
at a given pressure, temperature and velocity. Because of this, 
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FIG. 9: Size distribution of mesoscopic particle masses, measured 
in units of microscopic particle masses. The data was obtained from 
a simulation using 1600 microscopic and 160 mesoscopic particles, 
giving an average mesoscopic particle size of 10. 

the system is similar in spirit to the Lattice-Boltzmann coarse- 
graining, but with dynamic cells. An advantage of this method 
is that the dissipative part of the evolution equations can be de- 
rived theoretically 12411 . However, this method involves keep- 
ing track of, and updating, the Voronoi lattice at each time 
step of simulation, rendering this technique much slower than 
standard DPD. As the construction of the Voronoi lattice de- 
pends explicitly on all particle positions in the simulation, it 
also introduces a need for higher order interactions than the 
simple pairwise central forces normally associated with DPD. 
If it proves impossible to find a projection giving rise to DPD 
dynamics on the coarse grained level (which is a question that 
calls for further investigation), this alternative approach might 
still provide a reasonable path to take for performing reliable 
mesoscopic simulations. 



V. CHOICE OF PROJECTION 

Although the projection used in this study, Eqs. ( fT8l ), seems 
like a natural choice for DPD, it is a problem that the posi- 
tions of the coarse-grained particles only weakly reflect the 
positions of the underlying particles. As is shown in Fig. |9j 
the number of microscopic particles per coarse-grained parti- 
cle (i.e. the mass of the coarse-grained particle) exhibits large 
fluctuations, in sharp contrast to the standard DPD model 
where the masses of all particles are assumed to be equal and 
constant in time. 

One way to make the coarse-grained particles more closely 
reflect the density variations in the underlying system of mi- 
croscopic particles is to change the projection to incorpo- 
rate movement of coarse-grained particles towards regions of 
higher particle concentrations. This can be achieved by us- 
ing for instance the standard fc-means clustering method ll43ll 
(or any other position based clustering algorithm) to calcu- 
late the positions of the coarse-grained particles given the po- 
sitions of the underlying particles. This results in a model 
where the coarse-grained particles can be seen as clusters of 
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FIG. 10: The figure illustrates how the fc-means clustering algorithm 
works. A minimum in the potential landscape corresponding to a 
local minimum in the sum of squares function, Eq. i22i is found (in- 
dicated by vertical line). The small dots represent particle positions, 
Xi in Eq. i22\ . and V(c) is the sum of squares distance as a func- 
tion of cluster centre position. This example is made for illustrative 
purposes and therefore contains only a single cluster centre to which 
all the particles belong. As the particles move, V(c) changes, with 
the effect that local minima are continuously created and destroyed. 
This process results in discontinuous trajectories for the cluster cen- 
tre positions. 

underlying particles, with each cluster centre representing a 
local concentration peak of microscopic particles. An imple- 
mentation has been made using this projection, and the results 
reveal some new difficulties not easily foreseen in advance. 
With this type of projection, the cluster centres move in a po- 
tential landscape of the kind depicted in Fig. [10] for a one- 
dimensional system, where each local minimum of the curve 
represents a possible cluster centre position. The simulation 
was made for illustrative purposes, with a single cluster centre 
in a one-dimensional box with N = 100 particles, and with 
periodic boundary conditions. The curve represents the sum- 
of-squares distance from all the particles to the cluster centre, 
i.e. 

N 

V = J2 minflc - Xi\ 2 , \L-c+ x t \ 2 ) , (22) 

i=l 

where the minimum of the distance between the cluster and 
all periodically displaced images of particle i is used. 

By differentiating Eq. ( f22b with respect to the cluster cen- 
tre position, c, it is easily shown that a minimum in the sum 
of squares function represents a local average of the positions 
of the microscopic particles belonging to the cluster. This in- 
formation is just what the fc-mean clustering algorithm uses 
to calculate the cluster centre positions. The fact that the ex- 
ample in Fig. [10] contains only one cluster to which all the 
particles belong does not change the qualitative outcome that 
several local minima exist in the potential landscape. The re- 
sult of this is inevitably that the cluster centre positions, rep- 
resented by a given minimum in the potential landscape, will 
move with that local minimum until it disappears, which hap- 
pens frequently in the course of the simulation, for instance by 
the merging of two originally separated minima. At this point, 
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the cluster centre will jump to the adjacent minima, and in do- 
ing so it affects the neighboring cluster centres, resulting in 
discontinuous movement of the coarse-grained particles. As 
discontinuous particle movement on the coarse-grained level, 
due only to strictly local interactions on the microscopic level, 
is highly unsatisfactory, this type of position based projections 
also leave much to be desired. 



VI. SUMMARY 

In this article we have developed a method for estimating 
the forces between particles in a system that evolves accord- 
ing to the DPD ansatz, i.e. Langevinian dynamics with pair- 
wise central forces. The method works well for estimating 
both conservative and dissipative forces (with the stochastic 
force given by the dissipative through a fluctuation-dissipation 
theorem), and should work on any system that follows the 
DPD ansatz, as long as the time scale is small enough to let 
uj(r) be estimated from the force covariance k?. When ap- 
plied to a projected dynamics of a Lennard- Jones system, we 
cannot conclude that the projection results in a DPD-like dy- 
namics. The result points towards two possibilities: Either 
the projected dynamics is DPD-like, but outside the reach of 
our method, or in the worst case, there might be problems 
considering DPD as the result of a systematic coarse-graining 
method. 

A natural extension of the work presented in this article is to 
examine systems where artifacts due to fluctuating mass and 
identity problems are not encountered, such as the frequently 
used united atoms approach. A simple example would be to 
coarse-grain water by letting the coarse-grained particle be the 



whole water molecule. Some work in this direction has al- 
ready been made by Praprotnik et al. l44ll . Another direction 
is to develop a more sophisticated method for estimating cu(r) 
from Kp. 

As we suspect that the linear ansatz for is too simple, one 
might be tempted to simply use polynomials of higher degree 
in St and do a regression for the coefficent for each value of r 
separately, based on the region where we think the DPD the- 
ory is valid. Since Kp is close to linear in this region, however, 
the result of extrapolating the resulting function to find the in- 
tersection with the St = axis may be rather sensitive to the 
precise choice of region in St, and to noise in the measurement 
of Kp (from the finite number of samples). 

Rather, one may consider going in the other direction: for 
a given choice of ui(r) (and keeping the conservative force 
fixed) we measure ftp as a function of St and r and calcu- 
late a distance between ftp from the DPD simulation and Kp 
from the microscopic simulations. We may then use some op- 
timization procedure that does not require explicit calculation 
of derivatives, e.g. the classic Downhill Simplex method or 
Monte Carlo methods, to obtain better estimates for ui(r) (see 
e.g. i45ll for a review of different suitable optimization meth- 
ods). 
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